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ABSTRACT 

We present the first high resolution MHD simulation of cosmic-ray electron reacceleration by 
turbulence in cluster mergers. We use an idealised model for cluster mergers, combined with 
a numerical model for the injection, cooling and reacceleration of cosmic-ray electrons, to in- 
vestigate the evolution of cluster scale radio emission in these objects. In line with theoretical 
expectations, we for the first time, show in a simulation that reacceleration of CRe has the po- 
tential to reproduce key observables of radio halos. In particular, we show that clusters evolve 
being radio loud or radio quiet, depending on their evolutionary stage during the merger. We 
thus recover the observed transient nature of radio halos. In the simulation the diffuse emis- 
sion traces the complex interplay between spatial distribution of turbulence injected by the 
halo infall and the spatial distribution of the seed electrons to reaccelerate. During the forma- 
tion and evolution of the halo the synchrotron emission spectra show the observed variety: 
from power-laws with spectral index of 1 to 1.3 to curved and ultra-steep spectra with index 
> 1.5. 
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1 INTRODUCTION 

Galaxy clusters form in the knots of the cosmic web through the 
infall and merging of smaller structures. The gravitational poten- 
tial of these objects is dominated by dark matter (DM). The intra- 
cluster-medium (ICM) forms inside this potential well as a complex 
mixture of thermalised plasma and non-thermal components - mag- 
netic fields and cosmic-rays. Non-thermal emission is observed in 
the radio band from clusters (i.e. |Feretti et al.||2012"} , most spec- 
tacular in the form of giant radio halos. The origin of the underly- 
ing synchrotron bright cosmic-ray electrons (CRe) is still unclear, 
although the two commonly discussed mechanisms are reaccelera- 
tion by merger-driven tu rbulence l[Brunetti et al.||2001| |Petrosian 



2001 



2005 



|Ohno et al.||2002"l [FujiFa et al.||2003| |Cassano & Brunetti 



and the in-situ production of CRe by proton-proton collisions 



dDennison 



1980) |Blasi & C olafrancescol [T999l |Dolag & Englm] 



|2000||Pfrommer & En Blin 2004 1 or a combination of the two mech 
anisms ( Brunetti & Blasi 2005 ; Brunetti & Lazarian 201 1 1. To date, 
pure secondary models appear disfavoured by radio spectra of some 



halos ( Schlickeiser et al. 1 987 ; Thierbach et al. 2003 , Reimer et al. 
20031 |Brunetti et al.| '2008; Ma cario et~al 2010, Donner t~et~aT. 



2010 1. In addition 7-rays are unavoidably produced by the same 
decay chain as the CRe in these models, but to date, clusters re- 
main not observed in this regime {Aharonian 2009, Kiuch i et al.| 
|2009[ Ackermann et al. 2010}. This leaves hadronic models as- 



sailable, assuming cluster magnetic fields derived from the anal- 
ysis of Faraday rotation measures ( |Jeltema & Prof umo 201 1 ). Sev- 
eral works attempt to circumvent these difficulties by using a wider 
range of model parameters (Veritas Collaboration et al. 2012, Zan- 
|danel et al.[2 012). However a more detailed analysis that combines 
self-consistently state-of-the-art radio observations of the Coma ra- 
dio halo with the corresponding 7-ray (FERMI) upper limits chal- 
lenges a pure hadronic origin of CRe in this prototypical radio halo 
( [Brunetti etal.|2012) . 

Simulations of galaxy clusters have been used for a decade 
now to study the evolution of turbulence on the largest scales 
(Dolag et al. 2001 ). This has been done in the cosmological (Dolag 



et al.|2005[|Nagai et al.|2007l|Ryu et al.|2008|[Vazza et al.|2009 



,|2009| 



Iapichi no et al.||2011| l and"the idealised context~ ]Roettiger et al.| 
1999; Takizawa 2005). Here specifically idealised simulations of 
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cluster mergers usually neglect substructures of the two systems 
and the background expansion of the universe, but merge two 
spherically symmetric halos (Ricker & Sarazin 1 |200T] [Takizawa| 
|2008[ |ZuHone et al.|20lT) . These simulations allow for full con- 
trol of the system parameters and present a viable step towards 
fully cosmological approaches, which offer less control but more 
realism. 

Simulations of CRs in galaxy clusters have so far been limited 
to studies involing CR protons and the injection of CRe in shocks 
( |Miniati et al.|2001||Pfrommer et al.|2008l|Vazza et al.|201 1|[20T2> . 
These simulations did not include a treatment of turbulent reaccel- 
eration and were in fact unable to explain the main observed prop- 
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erties of giant radio halos simultaneously, mostly because of the 
underlying CR model i jDonnert et al. 2010 ?, B rown & Rudnickj 
|20TTT >. 

In this paper, we for the first time attempt to include the com- 
plex physics of CRe reacceleration by merger driven turbulence 
into simulations. This must be considered a first step in this direc- 
tion, as we model only CRe and simulate only an idealised binary 
collision. Still we are able to follow the evolution of turbulence, 
magnetic fields, CRe and the subsequent radio synchrotron emis- 
sion during mergers. We use a model for idealised cluster colli- 
sions based on an analytic solution of the hydrostatic equation in 
a Hernquist potential (Donnert in prep.). A low viscosity, high or- 
der description for SPH allows us to simulate the flow with high 
Reynolds numbers to follow the rise and decay of turbulence in the 
cluster. We then compute the evolution of CRe in post-processing 
to the simulation considering all relevant losses as well as for the 
first time stochastic reacceleration of CRe. Synthetic observations 
are then extracted using a numerical solver for the synchrotron in- 
tegral. 

This paper is organised as follows: in section [2] we summa- 
rize the CRe model assumed in our simulation. Section[3]describes 
the numerical method used to integrate the Fokker-Planck equation, 
treat turbulence in SPH and build the initial conditions. We present 
our results in section|4]and discuss these in|5] 



2 COSMIC-RAY MODEL 

In this paper, we consider only CRe and assume that they are re- 
accelerated by large scale MHD turbulence following |Brunetti "&] 
Lazarian ( 2007 ). There is agreement on the fact that turbulence in 
the ICM cannot efficiently accelerate thermal particles JPetrosian] 
|& East|2008) . However that mechanism is efficient enough to re- 
energize relativistic seed electrons (Petrosian 2001 ; Brunet ti~et al.| 
|200Tj|Ohno et al.|2002| >. Consequently, following previous works 
in the field, we assume a population of seed non-thermal electrons 
supplied by Q c (p,t). As this is an explorative study, we do not 
elaborate a physical model of the injection process of CRe. Possi- 
ble injection processes are diffusive shock acceleration, hadronic 
processes, reconnection or starbursting galaxies and AGN ( |Blasi| 
|et al.|2 007 ; Laza rian & Brunetti|201 1) . It has been realised early 
on that spatial diffusion of CRe is a slow process (see e.g. |Blasi| 
|et al.|2007| for a review). We therefore neglect it in this work. 

Considering an isotropic spectral energy distribution of CRe, 
n(p), the equation governing the evolution of CRe is then: 
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where 4| are the standard radiative and ionisation losses (Cassano 
|& Brunetti|2"005) : 
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where 7 = p e /m e c is the Lorentz factor and n t h the number den- 
sity of thermal gas. Specifically, we model the injection coefficient 
as: 

Q e ( 7j t) = ^ eet 2 h7 - 2 F (7) , (4) 



Value 


Unit 


Cluster 


Cluster 1 


P 




2/3 


2/3 


^Hcrnq 


kpc 


811 


473 


7-200 


kpc 


2192 


1380 


7*core 


kpc 


237 


130 


M2OO 


1O 15 M 


1.33 


0.17 


P0,gas 


10- 26 g/cm 3 


0.9 


1.1 




10 7 K 


9.2 


2.8 



Table 1. Model parameters of the initial conditions based on Donnert in 
prep. 

where eth is the thermal energy density of the plasma and F (7) 
is a cut-off function that limits the injection to energies of 7 6 
[50, 10 5 ] . This leaves K e the only free parameter in equation [T] 
Eq. [4] would eventually mimic continuous injection of secondary 
particles by inelastic collisions between CRp and thermal protons, 
with CRp following the spatial distribution of the thermal matter. 
For D pp = in equation[T] this results in the asymptotic solution: 

1 
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that is n(p) oc p~ 3 at high energies, where losses are dominated by 
synchrotron radiation and IC scattering (eq|2]l. To be conservative, 
we use only one mechanism to reaccelerate particles in the turbu- 
lent ICM, that is the interaction with fast modes due to the Transit- 
Time-Damping resonance (Schlickeiser & Miller] [1998). In this 
limit, the reacceleration coefficient considering the relevant colli- 
sionless damping of compressive turbulence in the ICM ( |Brunetti| 
& Lazarian 2007 1 can be simplified to : 
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where we assume a fraction 77 = 0.45 of turbulent energy in fast 
modes, c s is the sound speed, and «turb is me velocity of turbulenQ 
eddies at the scale I. 



3 NUMERICAL METHOD 



We use the MHDSPH code GADGET3 (Springel 2005) |Dolag &| 
|Stasysz yn 2009 1 to follow collisionless and gas dynamics of DM 
and the ICM plasma in the MHD approximation. 



3.1 Initial Conditions 

We model the collision of two dark matter halos using Hern- 
quist profiles, which can be identified with a NFW profile in the 
core ([Springel & Farrar 2007). The system has a total mass of 
1.5 x 10 lb M@, with a mass ratio of 1:8. We chose a compara- 
tively low mass-ratio in this simulation to reduce the influence of 
the shock phase to the radio emission. We leave a more general in- 
vestigation to future work. The density profiles are sampled by DM 
and SPH particles, 10 7 each. From the mass M200 in r2oo we obtain 
the Hernquist scale length using the observed relation of the con- 
centration parameter cnfw ~ 3 puffy et al.|2008) . For the bary- 
onic matter, we adopt the widely used /3-model, with f3 = 2/3 and 
fWe == 0.3r2oo/cNFW with a baryon fraction of 0.17 iRicker & 



|Sarazin 2001 , Mastropietro & Burkert 2008 1. For these parameters, 
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Figure 1. Velocity power-spectra of the inner most region of the system at 
different times. Top: after maximum halo brightness, bottom: before maxi- 
mum halo brightness. Time evolution is color coded. The dashed lines fol- 
low the Kolomogorov scaling of fc" 11 / 3 . 



the hydrostatic equation can then be solved analytically (Donnert 
in prep.). Similar to the observations of non cool-core clusters, the 
temperature profile is then almost constant in the central part and 
mildly declines towards the outer parts, with a mean temperature 
within the core region as listed in table[TJ alongside other details of 
the system. 

The two clusters are set-up separately, without magnetic field, 
and relaxed for several Gyrs. The relaxation process is done with 
high viscosity and thermal conduction (Jubelgas et al. 2004|>, to 
quickly equalise density fluctuations induced by the Poisson sam- 
pling of the gas density distribution. 

In a second step, the clusters are joined in a periodic box of 
10 Mpc. This allows a consistent treatment of the magnetic field 
in k-space. Furthermore, periodicity prevents low timesteps on es- 
caping SPH particles. The two halos are put on a zero energy or- 
bit at a distance so that the density of both profiles match with 
an impact parameter of 300 kpc. Overlapping particles are resam- 
pled to the volume outside of V2oo- The magnetic field is then ini- 
tialised as a constant divergence free vector-field in k-space, with 
a Kolmogorov power-spectrurr[^]down to 150 kpc. This minimises 
the oversampling of the SPH kernel in low density regions, which 
is known to cause numerical divergence. The vector-field is then 
transformed to real-space and sampled to the particles with the TSC 
kernel. The resulting field is rescaled to a maximum value of 3 /j,G 
in the cluster centre and attenuated radially as the square-root of 
the gas densit)|^] Residual magnetic divergence is then treated by 
the divergence cleaning in the code at runtime (see section[3~2]l. 



3.2 Treatment of Turbulence 

To follow the turbulent reacceleration of the CRe, we need to trim 
the SPH formalism to follow the generation of turbulent motions 
and their propagation, possibly over a broad range of scales within 
the simulation. Furthermore, we have to estimate the turbulent en- 
ergy at every location within the cluster during the merger event. 
Therefore, we use a largely improved SPH scheme, where: 



Kuchar & EnBlin 



2011 



as measured by e.g. 
as derived from Faraday rotation measurements by Bonafede et al. 2010 



in Hy dra A 



(i) Instabilities can be resolved by using a time-dependent for- 
mulation of thermal conduction (Price 2008), which allows a mix- 
ing similar to grid codes. This additionally eliminates discontinu- 
ities in our initial conditions. 

(ii) We employ a time-dependent viscosity scheme based on the 
local signal velocity (e.g. |Dolag et al.|2005) , which significantly 
reduces the numerical viscosity outside of shocks. 

(iii) We use a higher order kernel based on the bias-corrected 
fourth order Wendland kernel pehnen & Aly||2012) with 200 
neighbors. This corresponds to the same effective kernel FWHM 
(smoothness) as the standard, cubic kernel with 64 neighbours, but 
increases the compact support of the kernel considerably (Dehnen 
|& Aly|2012^ . Then dissipation takes place almost entirely inside 
this kernel scale. 

This way the SPH algorithm formally models a flow, with 
a Reynolds number of more than 1000 away from shocks, and 
sufficiently follows turbulent motions up to the kernel scale. In 
such an implementation, SPH does form a cascade down to small 
scales (Hopkins 2012), where turbulent motions are eventually 
thermalised and damped by the viscosity scheme. This approach 
is therefore conservative in terms of turbulent energy. 

In figure[TJ we show the filtered velocity power- spectrum of a 
region of 1 Mpc 3 around the center of mass of the system before, 
during and after the merger. Here we binned the SPH particle ve- 
locity to a grid. Then we subtracted the mean velocity in a region of 
100 kpc around each cell to filter out the bulk motion of the system. 

The intial power-spectrum shows fluctuations of the order of 
lOkm/s, due to the underlying numerical algorithm. The amplifi- 
cation of the velocity power during the passage is indicative of a 
drastic increase in turbulent motions in that phase. Furthermore, the 
subsequent decay of that velocity power then demonstrates the tur- 
bulent origin of these motions and their dissipation on small scales. 
Thus, we use a high order kernel which employs the same effec- 
tive size, but a much larger region and number of neighbors than 
standard approaches. This allows us to estimate the local turbulent 
energy, similar to sub-grid turbulence models for modelling unre- 
solved turbulence in galaxy clusters llapic hino et al.p Olll. Here 
we use the RMS velocity fluctuation filtered on the kernel scales. 
This implies I — 2/i sm i in equation [6] where h sm \ is the compact 
support of the SPH kernel. Note that compared to grid codes, this 
method typically underestimates the turbulent energy on the kernel 
scale. 

We model a turbulent magnetic resistivity, following Bonafede 
|et al.| ( |2011| l. This also acts as a cleaning scheme for magnetic di- 
vergence in turbulent flows and treats residual divergence in the 
outer regions as well as the initial conditions (Pric e|20T2) . 



3.3 Cosmic-Ray Electrons 

We solve equation [T] for every SPH particle, in post-processing to 
the simulation. We employ the adaptive upwind scheme developed 
by Chang & Cooper ( 1970). Per particle, we logarithmically sam- 
ple the isotropic number density spectrum of CRe 100 times in the 
range of p/m e c g [0.1, 10 5 ], with open boundary conditions re- 
alised, as in Borovsky & Eilek ( 1986). The fluid quantities are in- 
terpolated linearly between 290 snapshots every 25 Myr. We em- 
ploy a universal timestep of 0. 1 Myr. For example, we show the 
spectrum of one particle in figure|2] It shows the effect of turbulent 
acceleration on the spectrum of the CRe. Initially (dark blue), the 
system shows the asymptotic solution (eq.|5j, then (light red, red) 
turbulent acceleration modifies the spectrum up to 7 ^ 3 x 10 4 . At 
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Figure 2. Spectrum of an exemplary particle (ID 400383) at different times 
during the merger. Colours correspond to figure[4]and[5] 



the highest energies, the spectrum is shaped by the constant injec- 
tion of particles (eq.|4|. The norm of this injection varies with time 
as the underlying thermal conditions evolve. At late stages (green, 
light green, light blue), turbulent acceleration is less efficient, lead- 
ing to the accumulation of particles at synchrotron dark energies 
around 7 < 10 3 . 

To produce synthetic observations from our simulations, we 
use our MPI-parallel projection code SMAC2, and numerically 
solve the synchrotron integral for every particle. 



4 RESULTS 

The simulation is evolved for roughly 4 Gyr, and CRe spectra and 
synthetic observations are obtained. The first core-passage occurs 
at ~ 1.18 Gyrs. In figure[3] we show X-ray emission of the system 
at different times. 

Despite the simulation being highly idealized, which is re- 
flected in morphological details (high symmetry and regular shapes 
in the X-rays), the system can be classified into relaxed (before) 
and disturbed (after) state. We find that Mpc-scale synchrotron ra- 
dio emission is generated as a result of particle acceleration by 
merger turbulence, in connection with the collision between clus- 
ters. The radio - X-ray evolution of the system can be summarized 
in 3 phases: infall and shock dominance, reacceleration phase and 
decay. 

(i) During infall (less than 1.25 Gyr after the start of the simu- 
lation), X-ray and radio luminosity increase, within a few hundred 
Myr, by a factor of 2 and 30, respectively. The radio emission is 
localized at the region crossed by the shock front. 

(ii) The system enters the reacceleration phase shortly after the 
core passage (1.25 < t < 1.7 Gyr): The X-ray brightness rapidly 
declines due to the decrease in density, caused by the disturbance of 
the core. In contrast, the total radio synchrotron brightness contin- 
ues to increase, as the DM core drives turbulence in a large fraction 
of the cluster volume. 

(iii) During the decay (t > 1.7 Gyr) phase, the radio emission 
traces the trail of the turbulence driven by the secondary DM core, 
fades away gradually and becomes offset from the primary core. 

Present radio surveys suggest that radio halos are transient 
phenomena connected to cluster mergers (Cassano et al. 20l0|>. In 
figure]?] we plot the evolution of the system in the Lx - P14 plane as 



Figure 4. Evolution of the system in X-rays and radio synchrotron emission 
as black line. We mark the time evolution as dots on the lines at an inter- 
val of 25 Myr. The colored dots correspond to the times shown in figure 
[3] The dashed line shows the emission from the injection only. Overplot- 
ted are observed radio halos and upper limits tBrunetti et al. 2009} as well 
as the correlation (black line). To guide the eye the correlation is scaled to 
Coma (black dotted). We mark Coma with a black cross. We also include 
known upper limits as arrows. Ultra steep spectrum halos A521, A1914, 
MACSJ1 149.5+2223 and A697 are plotted as asterisks (Dallacasa eTaf] 
pO^[vanWeeren et al.|2"0TT||Bonafede et al.|2012) . 



a black line. This is compared with observed radio halos (astrisk - 
ultra-steep spectrum) and present upper limits. We mark time inter- 
vals of 25 Myr as dots on the line, and the times shown in figure[5]as 
large dots with corresponding colour. The emission from the injec- 
tion only, i.e. without considering turbulent reacceleration, is over- 
plotted as dashed line. To determine the CRe normalisationi, a time 
was chosen, so that the simulated spectrum fits the observed Coma 
spectrum. Then the normalisation was set to K a — 2 x 10 -4 , so that 
the halo brightness roughly fits the Coma luminosity extrapolated 
along the correlation (dotted line in figure [4j. This corresponds to 
a total energy injected in the form CRe during our simulation of: 



£CRe,inj 
Eth 
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2 Gyr 



< eth > 
2 x IO" 11 

cm° 



(7) 



is the 



where < eth > is the typical thermal energy density and n n] 
injection time. 

Before the first passage, the radio luminosity of the system is 
more than a factor of 10 below the correlation. This corresponds to 
the "off-state" of galaxy clusters without an observed radio halo, 
i.e. the upper limits from observations. During the shock phase of 
the first passage, X-ray and radio luminosity of the system increase 
and the system approaches the zone of radio halos. In the reaccel- 
eration phase, the system crosses this zone as the X-ray luminosity 
declines rapidly by a factor of two. The radio luminosity decreases, 
within 1 Gyr, to injection values, i.e. the system enters the radio- 
quiet phase. The system spends ~ 0.7 Gyrs inside the observed 
scatter around the correlation and ~ 0.17 Gyrs to move between 
correlation and upper limits. 

The spectral properties of the system change drastically dur- 
ing its evolution in the P14-LX plane. In figure|5] we plot the ra- 
dio spectrum of the system at different times. We overplot the ra- 
dio spectrum of the COMA halo (see Thierbach et al. 2003] and 
ref. therein) as well as A521 (jDallacasa et al. 2009). The simulated 
spectra of Coma and A521 are scaled to match the observed ones 
at 430 MHz and 240 MHz, respectively. 

Before the merger, the system shows only emission with a 
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Figure 3. Projections of bolometric X-ray luminosity of the system rotated by 45. Overlayed radio contours (1.4GHz) are at 0.050, 0.22, 7.3, 1.45 X 10' 
erg s — 1 . 
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Figure 5. Radio synchrotron spectrum of the system at different times. Ob- 



served spectrum of the Coma cluster as black diamonds from |Thierbach| 
et al.|{2003 1. Observed spectrum from A521 as asterisks (Dallac asa et al.| 



20091. 



spectral index of ~ —1, consistent with the steady-state spectrum 
due to losses and constant injection (eq. [2] [3] |4j- At maximum lu- 
minosity, the radio spectrum fits the observed COMA spectrum. 
A comparison with A521 yields a reasonable fit of the spectrum 
at 1.8 Gyr after start of the simulation. Over time, the total emis- 
sion decreases and further steepens at high frequencies. Different 



assumptions for the fraction of compressive modes would lead to 
systematically flatter (high rf) or steeper (low 77) spectra. However, 
for the first time, we have shown that viable assumptions for this 
parameter in simulations enable us to naturally explain the variety 
of observed spectra ( Venturi 201 1 ) in giant radio halos. 



5 CONCLUSIONS 

We obtained, for the first time, an MHD simulation of two collid- 
ing, idealized galaxy clusters, including the evolution of the mag- 
netic field as well as a combined treatment of turbulence and CR 
electrons reacceleration within the Reacceleration is treated 

by a simple conservative model, assuming TTD damping of fast 
magneto-sonic waves by relativistic particles. The simulated evo- 
lution of the ICM and all relevant cooling mechanisms are cou- 
pled self-consistently to the CRe. We numerically solve the under- 
lying Fokker-Planck equation for all 10 7 gas particles of the cluster 
merger simulation, sampling the CRe distribution function by 100 
logarithmically-spaced energy bins, in the range 7 = [0.1, 10 s ] . 
This allows us to obtain synthetic radio observations in the ob- 
served frequency range (30MHz - 3GHz) and compare to the dis- 
tribution and evolution of the X-ray emission of the system. The 



4 We note that a similar approach was adopted for mini-halos and core- 
sloshing by ZuHo ntTet al.||2012| 
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inclusion of these techniques comprises a significant step towards a 
self-consistent treatment of CR physics in large-scale astrophysical 
simulations. 

Despite the idealised set-up of the underlying thermal model, 
our simulation for the first time reproduces several key observables 
of giant radio halos. In particular: 

(i) The transient nature of radio halos and their connection to 
mergers and merger driven turbulence. 

(ii) The variety of observed radio spectra (i.e. flatter, curved and 
ultra-steep spectrum halos) associated with different states in the 
merger evolution. 

Our results differ from those of previous simplified numerical mod- 
els, that were based solely on hadronic collision in the ICM and 
shock injection of CRe and CRp. 

These results reaffirm previous theoretical arguments that 
curved (or very steep) spectra of radio halos, their transient na- 
ture and connection with mergers are solid predictions of turbulent 
reacceleration models (Brunetti & Lazarian 2007, 2011). Present 
results must be considered to be only first steps in this direction, 
as a more self-consistent model requires the inclusion of CRp, the 
generation of secondary particles and the interplay with turbulence 
in a cosmological setup. 
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